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We propose a methodology to analyze synchronization in an ensemble of diffusively coupled multistable 
systems. First, we study how two bidirectionally coupled multistable oscillators synchronize and demonstrate 
the high complexity of the basins of attraction of coexisting synchronous states. Then, we propose the use of the 
Master Stability Function (MSF) for multistable systems to describe synchronizability, even during intermittent 
behaviour, of a network of multistable oscillators, regardless of both the number of coupled oscillators and the 
interaction structure. In particular, we show that a network of multistable elements is synchronizable for a given 
range of topology spectra and coupling strengths, irrespective of specific attractor dynamics to which different 
oscillators are locked, and even in the presence of intermittency. Finally, we experimentally demonstrate the 
feasibility and robustness of the MSF approach with a network of multistable electronic circuits. 

PACS numbers: 89.75.Hc,89.75.Fb 


INTRODUCTION 


During the last two decades great progress has been made in 
understanding synchronization of chaotic systems Hi . Nev¬ 
ertheless, synchrony in dissipative dynamical systems with 
coexisting attractors remains relatively unexplored and poorly 
understood to this very day. This relative lack of activity is 
hard to reconcile with the fact that multistability has been ob¬ 
served in numerous nonlinear systems in many fields of sci¬ 
ence, such as laser physics 0], neuroscience flfe], cardiac 
dynamics |@], genetics US], cell signaling fsjf], and ecology 
110] amongst others; moreover, in situations where synchro¬ 
nization is actually a collective behavior known to play a pri¬ 
mary role. Even forms of extreme multistability, i.e. the coex¬ 
istence of infinitely many attractors in phase space, have been 
recently observed in experiments 0. Many of these results 
as well as some known coupling mechanisms and dynamical 
phenomena that seem to be correlated to the emergence of 
multistability are reviewed in Il2il . 

The dynamics of two unidirectionally coup led systems 
(master-slave configuration) of Rossler-like ll 13[ 114fl, Duffing 
[15] and Rossler-Lorenz II1611 oscillators, as well as Henon 
maps 10 have been studied, and some experiments along 
these lines have been carried out 1811 . Bidirectionally cou¬ 
pled neuronal models 0 also display very rich synchronous 
dynamics. One of the most prominent features of all these 
examples is the intricate dependence of synchronization on 
the initial conditions, a distinct feature of multistable sys¬ 
tems that is nowhere to be found in monostable systems. 
Phenomena such as anticipated intermittent phase synchro¬ 


nization, period-doubling synchronization, and intermittent 
switches between coexisting type-I and on-off intermittencies 
have been discovered 1141 . 118 1. On the other hand, even though 
the existence and stability of multistable synchronous solu¬ 
tions in locally coupled Kuramoto models have been studied 
100, to our knowledge, the issue of under what conditions 
synchronization of more than two coupled generic (and possi¬ 
bly chaotic) multistable oscillators is guaranteed has not been 
addressed yet. Furthermore, synchronization of multistable 
systems in the presence of intermittency still remains an un¬ 
explored problem. 


In this paper, we propose a methodology for studying the 
synchronization of multistable oscillators, which we illustrate 
with the example of a bistable system which has the great ad¬ 
vantage of being experimentally implemented in electronic 
circuits. First, we demonstrate the high complexity of the 
basins of attraction of coexisting states in a solitary bistable 
oscillator, and the increasing complexity when two of such os¬ 
cillators interact with each other giving rise to intermittency. 
Second, we investigate the influence of both the initial condi¬ 
tions and the coupling strength on the synchronization of two 
bidirectionally coupled bistable systems (with diffusive cou¬ 
pling) in different coexisting synchronous states, including the 
existence of intermittency. Then, we discuss the Master Sta¬ 
bility Function (MSF) |23[] approach to the study of the sta¬ 
bility of a synchronization manifold of N coupled multistable 
systems. Specifically, we obtain the MSF for different coex¬ 
isting chaotic attractors in a dynamical system separately, and 
then we evaluate how the modification of the coupling param¬ 
eter allows the system to leave/enter a particular synchroniza- 
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tion regime associated with a particular attractor without loss 
of synchrony in the whole network, even in the presence of 
intermittency. Finally, we check the robustness of our theo¬ 
retical predictions with electronic circuits to show the validity 
of our results for real systems where a certain parameter mis¬ 
match always exists. 


CHARACTERIZATION OF THE SYSTEM IN TERMS OF 
ATTRACTORS AND THEIR BASINS OF ATTRACTION 


The main aim of our work is to develop a methodol¬ 
ogy that adequately predicts synchronizability in ensembles 
of bidirectionally coupled multistable systems. With this 
aim, we choose the piecewise linear Rossler oscillator as a 
paradigmatic example of a bistable system with two coexist¬ 
ing chaotic attractors |Q. When an unidirectional coupling 
is introduced, the coupled Rossler-like system exhibits very 
rich dynamics, including such phenomena as intermittency, 
frequency-shifting and frequency-locking [13;]. Nevertheless, 
little is known about synchronization scenarios for ensembles 
of bidirectionally coupled multistable systems, despite the fact 
that bidirectional coupling itself can lead to the emergence of 
multiple attractors [24]. 

Specifically, the equations describing the dynamics of the 
Rossler-like oscillators are i ll311 


x = —a\(x + /3 y + I\s), 
y = -a 2 (-'yx - [1 - <%), 
z = - a 3 (-g(x ) + z), 


(1) 


with 


/ \_ J 0 x < 3, 

i n (x-3) x > 3, 


( 2 ) 


where x, y, and z are the state variables. The piecewise lin¬ 
ear function g(x) introduces the nonlinearity in the system 
that leads to a chaotic behavior. The parameter values are 
Qi = 500, a 2 = 200, c*3 = 10000, /? = 10, T = 20, 7 = 50, 
5 = 15.625 and y = 15. For this parameter choice, the system 
is known to be a bistable chaotic system (i.e. the phase por¬ 
trait displays two different chaotic attractors), as previously 
reported in d- Unlike most previously studied multistable 
systems (see, e.g. II16|] ). this system exhibits multistability in 
the autonomous evolution, without the need for chaotic driv¬ 
ing. Moreover, this system can be implemented in electronic 
circuits to experimentally assess the validity of the theoretical 
predictions. 

From any arbitrary initial condition within a bounded re¬ 
gion in phase space the system rapidly converges to one of 
the two chaotic attractors shown in Fig. (T|(a). We denote the 
larger attractor by L (Fig. [j] (a), blue (dark gray)) and the 
smaller one by S (Fig. |T|(a), red (gray)). The basins of at¬ 
traction of L (blue (dark gray)) and S (red (gray)) are shown 
in Fig. ffl(b) for initial conditions x(0) = (x(0), y(0), z( 0)) 
such that x(0) £ [—4,6], y( 0) £ [—8,4] and ,z(0) = 0. The 



FIG. 1 . (Color online) Coexisting attractors and their basins of 
attraction in the Rossler-like oscillator (Eq. |T}. (a) Large ( L, blue 
(dark gray)) and small (5, red (gray)) 
attractors, (b) Basins of attraction for L (blue) and S (red) in the 
z = 0 plane. Initial conditions leading to unstable trajectories 
appear in white, (c) Basins of attraction for L and S, 

[—1,1] x [—1,1] square in the z — 0 plane, (d) Basins of attraction 
for L and S, [—0.1, 0.1] x [—0.1, 0.1] square in the z = 0 plane. 


basin of attraction of L is seen to be much larger than the basin 
of S. Two spirals are clearly visible, where initial conditions 
leading to one or the other attractor seem to be intertwined. 
Each of these spirals has a fixed point of the system as its fo¬ 
cus, as has been previously reported in 0 . In Fig. 00(c) 
we focus on x( 0 ) £ [— 1 , 1 ], y{ 0 ) £ [— 1 , 1 ] to better appre¬ 
ciate the details of the spiral that has its center in the origin. 
Indeed, the mixing of initial conditions close to the center of 
the spiral seems to be present at arbitrarily low space scales, 
as the zoom around x( 0 ) £ [— 0 . 1 , 0 . 1 ], y{ 0 ) £ [— 0 . 1 , 0 . 1 ] in 
Fig. 00(d) shows. We have checked that this structure is pre¬ 
served for another 4 orders of magnitude, with no end in sight 
for even lower space scales. Indeed, this is not altogether sur¬ 
prising, as basins that are interwoven in a complicated fashion 
and fractal basins boundaries feature prominently in the phase 
portraits of many multistable systems (see li25N for review of 
fractal basin boundaries and fractal sets in nonlinear dynam¬ 
ics in general). Although the precise characterization of these 
boundaries is outside the scope of this paper, we would like to 
stress how difficult it is to control the asymptotic dynamics of 
just one single oscillator under the presence of noise or uncer¬ 
tainties for initial conditions starting in certain regions of the 
phase space. 


TWO BIDIRECTIONALLY COUPLED SYSTEMS: 
COEXISTING REGIMES 


To start our study on synchronization of bidirectionally 
coupled multistable systems, we first consider the simple case 
of two coupled Rossler-like oscillators. This particular sys- 
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tem has been thoroughly analyzed for the case of a master- 
slave configuration ll3l 114111811 . Here, we investigate, for 
the first time to our knowledge, bidirectionally coupled multi¬ 
stable chaotic systems, and also for the first time the coupling 
is diffusive. The coupling is introduced through the x variable 
with coupling strength a so that the equations of the motion 
become 

it, 2 = -ai(ati ,2 + Pyi ,2 + I>L ,2 - CTIp (£ 2,1 - £1,2)), 
2/1,2 = -ct 2 (—7 £i,2 - [1 - <%1,2), 
il,2 = -« 3 (-ff(£l, 2 ) + 21,2), 

(3) 

where ip = 20 and g(x) is given by Eq. (0. Our numeri¬ 
cal simulations show the existence of four possible asymptotic 
regimes: a) both systems end up in an attractor indistinguish¬ 
able from L (from now on, LL), b) both systems end up in an 
attractor indistinguishable from S (SS), c) one system asymp¬ 
totes to L and the other to S (SL ), d) the systems switch inter¬ 
mittently back and forth between L and S in an irregular way 
(!) (the intermittent behavior of these systems in master-slave 
configurations is described in [11811'). 

As the coupling strength a is increased from 0 to 0.2, all 
these cases appear, disappear and mix in a very complicated 
manner, depending on the initial conditions. In Fig. [2] we 
show the basins of attraction for these asymptotic regimes that 
result from fixing 2 /i( 0 ) = 2 / 2 ( 0 ) = 0 and 2 i( 0 ) = 22 ( 0 ) = 0 , 
and exploring a finely discretized grid for £ 1 , 2 ( 0 ) £ [—4,4] 
for four representative coupling strengths, er = 0, 0.05, 0.14, 
and 0.2. For moderately larger er the basins remain relatively 
stable, while the phase space projections of the full attractor 
on the subspace corresponding to each of the two subsystems 
display only very slight deformations with respect to the orig¬ 
inal L and S seen for er = 0. 

For better visualization of the prevalence of the different 
asymptotic regimes for different values of er. Fig. [3] (a) shows 


(a) (b) 



FIG. 3. (Color online) Relative size of the basins for each attrac¬ 
tor, synchronization error, largest and second largest Lyapunov 
exponents of the full system as a function of a. (a) Relative size of 
the basins of attraction of the full system, (b) synchronization error 
(inset: zoomed part of the plot), (c) largest Lyapunov exponent, (d) 
second largest Lyapunov exponent. 


the fraction of the initial conditions considered, which even¬ 
tually lead to every possible asymptotic regime, as a function 
of the coupling strength. As expected from Fig. □ Fig. 0(a) 
shows that for small a most of the initial conditions lead to the 
intermittent behavior (Fig. 0(b)), whereas stronger coupling 
leads first to the emergence of SS (Fig. 0 (c)) and later to 
the emergence of LL together with a vanishing of intermittent 
behavior I. For large cr, only LL and SS are observed (Fig. 0 
(d)). Interestingly, the SL case is never observed for nonzero 
a, so that we do not consider it anymore. 



FIG. 2. (Color online) Basins of attraction of the system compris¬ 
ing two bidirectionally coupled Rossler-like oscillators for differ¬ 
ent coupling strengths a. The coupling strengths are (a) a = 0 
(uncoupled systems), (b) a = 0.06, (c) er = 0.14, and (d) a = 0.20. 


TWO BIDIRECTIONALLY COUPLED SYSTEMS: 
SYNCHRONIZATION OF DIFFERENT REGIMES 


Having seen the different possible asymptotic regimes and 
their basins of attraction, we move on to study synchroniza¬ 
tion of two bidirectionally coupled systems. Instead of focus¬ 
ing on individual trajectories, we average the results across all 
initial conditions leading to one of the considered cases (LL, 
SS or I). First, we calculate the synchronization error (Fig. 
0(b)) for each of the different asymptotic regimes. This way, 
we observe how synchronization is first achieved for the SS, 
red (gray) triangles, case around a ~ 0.075, and for the LL, 
black (dark gray) diamonds, case only for a ~ 0.165. As 
for the I, green (light gray) squares, case, the synchronization 
depends on the coupling in a more complex, non-monotonic 
manner, which in some way reflects the synchronization pat¬ 
terns of SS and LL. As the system becomes more and more 
synchronized and therefore the influence of one subsystem on 
the other becomes weaker, I fades away. The transition to syn¬ 
chronization in this regime is consistent with the observation 
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made in 11511 . as it passes through two stages: first, the two 
oscillators move to the same attractor due to the increasing 
coupling, and second, as the coupling further increases, they 
gradually synchronize (just as monostable coupled oscillators 
do). 

A different view on the emergence of synchronization is 
provided by the analysis of the Lyapunov spectrum. Figured 
(c) shows the largest Lyapunov exponent as a function of a 
for different asymptotic regimes. The dotted lines show the 
largest Lyapunov exponents corresponding to the attractors L 
and S for the solitary (uncoupled) system. In the I regime, the 
value of the Lyapunov exponent is clearly correlated with the 
presence of LL and SS windows, a phenomenon that will be 
elucidated below. However, from the synchronization point of 
view, the most important information derived from the Lya¬ 
punov spectrum concerns the second largest Lyapunov ex¬ 
ponent, which gives indirect information on the existence of 
weakest forms of synchronization, such as generalized syn- 
chronization (GS) (see, for instance, the discussion in [26il). 
As seen in Fig. 0d), two points, highlighted as ass and gll, 
corresponding to the loss of the positive Lyapunov exponent in 
the system, mark the starting point of the existence of just one 
chaotic mode in the system, and therefore the onset of GS for 
SS and LL respectively. Also, I seems to reach GS, but as we 
will see I is a very sophisticated regime, with frequent jumps 
between time windows, where both oscillators are in the same 
attractor, and time windows, where each oscillator is in a dif¬ 
ferent attractor. Time averages across very long time windows 
(as those used in the calculation of the synchronization error 
or Lyapunov exponents) fail to capture this sort of details, so 
we presently conclude this section with a brief study of this 
regime to illustrate just how complex it can be. 

We now focus on the dynamics of the coupled system in 
the I regime and measure the fraction of time spent by both 
subsystems simultaneously in the S state (denoted as ss ), L 
state ( ll ), and when each subsystem follows a different dy¬ 
namics (si). Figure 0] (a) shows the fraction of time that the 
system being in the I regime spends in ss, red (gray) trian¬ 
gles, ll, green (light gray) squares, and si, black (dark gray) 
diamonds, time windows, which bears an obvious relation to 
the results shown in Fig. 0(c). Another important issue is the 
lack of complete synchronization in I even for very high a. 


(a) (b) 




c a 


FIG. 4. (Color online) Fraction of time and synchronization error 
in the / regime for different time windows as a function of a. (a) 

Fraction of the time spent by the system in ss, ll, and si windows, 
(b) I synchronization error restricted to ll and ss time windows. 


To examine this problem in more detail, we now compute the 
synchronization error separately for the ss and 11 time win¬ 
dows. We do not consider the si case because by definition 
complete synchronization is not possible in these time win¬ 
dows. In Fig. 0(b) the synchronization error reveals that for 
a > 0.075 the I regime contains completely synchronized ss 
time windows, while the ll time windows show a clear lack of 
synchronization. It is difficult to ascertain whether the mild 
growth in the synchronization error of ss for a > 0.10 is in¬ 
deed a loss of synchronization or has to do with the difficulties 
associated with the classification of time windows into ss, ll, 
and si & In the next section we will present a straightfor¬ 
ward way to understand the conditions under which synchro¬ 
nization is possible in an ensemble of multistable dynamical 
systems even in the presence of intermittent regimes such as 
the one we just described. 


MASTER STABILITY FUNCTION OF MULTISTABLE 
DYNAMICAL SYSTEMS 


Up to now, we have seen how complex the basins of at¬ 
traction of the bistable oscillators are, and how they become 
much more complex as soon as two oscillators interact, with 
the development of new phenomena such as intermittency, 
that makes synchronization predictions very difficult and ex¬ 
tremely dependent on initial conditions. Whenever multi¬ 
stable oscillators are coupled in a large ensemble, we expect 
the dynamics to become even more complex, and the prob¬ 
lems associated with the study and control of synchronization 
become even more serious. In order to measure the synchro¬ 
nization stability of ensembles of coupled multistable oscil¬ 
lators, we resort to the Master Stability Function (MSF) ap¬ 
proach of Pecora and Carroll |23tl . The MSF curve is given 
by the maximum of the Lyapunov exponents transverse to the 
synchronization manifold. In the past 15 years, the MSF ap¬ 
proach has found important applications in the study of net¬ 
works of monostable systems (see, e.g. |28|), however, as far 
as we know, it has never been extended to the case of systems 
with coexisting attractors. 

While the MSF describes the linear stability of the syn¬ 
chronous motion for a given attractor dynamics, in the pres¬ 
ence of multistability one should rather look at the MSF of 
each attractor separately. In Fig. 0(a) we show the MSFs cor¬ 
responding to both S, red (gray) continuous line, and L, blue 
(dark gray) dashed line, for the Rossler-like system with dif¬ 
fusive coupling through the x variable as in Eq. (0. The 
independent variable v is a parameter that implicitly takes 
into account infinitely many network topologies and coupling 
strengths, as it stands for the product of a given coupling 
strength a and a graph Laplacian eigenvalue (for each eigen¬ 
mode transversal to the synchronization manifold there is one 
such eigenvalue, and there are N — 1 of them for a fully con¬ 
nected network of size N). According to the classification 
proposed in 1281. this system belongs to the class III systems, 
as the region of stable synchronization is bounded between 
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two MSF zeros. Interestingly, Fig. [5] also shows that the sta¬ 
bility region of L is contained in that of S. The interplay be¬ 
tween the two stability regions ( S and L) given by the MSFs 
explains how the synchronized state is maintained (or lost) in 
one of the three possible scenarios SS, LL, and I. There are 
therefore three possibilities for a given a and a given eigen¬ 
mode of the topology (i.e. for a given v = crA, where A is 
the specific eigenvalue of the network’s Laplacian matrix that 
corresponds to the considered eigenmode): a) neither S nor L 
can be synchronized (i.e., the value of the MSFs of both S and 
L are positive for such ;/), b) S can be synchronized, but not 
L (i.e., the MSF is negative for S, but positive for L), and c) 
both S and L can be synchronized (both values of the MSFs 
are negative). While b) only guarantees the synchronization 
stability in the SS regime, c) is tantamount to saying that syn¬ 
chronization is stable no matter how complex the dynamics 
may be, even in the presence of intermittency. Therefore, if 
the product of each of the eigenmodes for a given topology 
and a is within the region described by c), synchronization 
is stable. Analogous arguments can be used in systems with 
an arbitrarily large number of attractors, provided the stability 
regions of the different attractors are not disjoint. 

Figure 0(b) shows the MSFs corresponding to the two co¬ 
existing attractors when the oscillators are coupled through 
the y variable. In this case the system belongs to class II, 
i.e. the stability region starts at a given v and then extends 
indefinitely to the right. Again, for this particular system the 
stability region of L is contained in that of S. Choosing a 
high enough (or changing suitably the topology for a given 
a) in such a way that all eigenmodes are in the region lying 
to the right of the MSF zero for L guarantees the stability of 
synchronous dynamics on any attractor, even in the presence 
of intermittency. Below we will dwell upon the contents of 
the insets in Figs. |5](a) and (b). 
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FIG. 5. (Color online) MSFs for the coexisting attractors for two 
types of coupling, and variability with respect to uncertainties in 
the parameters, (a) MSF of S, red (gray) continuous line, and L, 
blue (dark gray) dashed line, for oscillators coupled through the x 
variable (class III system); inset: variability of the first zero crossing 
as the parameters are affected by uncertainties (see Section Exper¬ 
imental Implementation for a detailed explanation), (b) MSF of S 
(red continuous line) and L (blue dashed line) for oscillators coupled 
through the y variable (class II system); inset: variability of the zero 
crossing as the parameters are affected by uncertainties. 



EXPERIMENTAL IMPLEMENTATION 


The MSF approach determines the stability of the synchro¬ 
nization manifold but does not give information about how 
large the stability basin of the system is. Nevertheless, in 
real systems, the size of the stability basin of the synchro¬ 
nized manifold plays a fundamental role in the observation 
of synchronization, specially when external perturbations are 
applied (see |[29] for a detailed description of how to evaluate 
a stability basin). In this section, we compare our theoreti¬ 
cal predictions given by applying the MSF approach to multi¬ 
stable systems with experimental results in order illustrate the 
robustness under the presence of noise and parameters mis¬ 
matches, and verify their validity even when the noise and the 
coupling of non-identical systems bring the system quite far 
away from the synchronization manifold. 

The experimental design is based on a 6-node network of 
piecewise Rossler-like electronic circuits coupled according 
to the topology shown in Fig. [6j he. following a spiderweb 
network topology with a central node connected to all other 


nodes, and each of the 5 peripheral nodes connected to their 
2 spatial neighbors. Other network configurations, with an ar¬ 
bitrary number of oscillators, are also possible, the only lim¬ 
itation being the number of available electronic components. 
The circuits are the experimental implementation of the dy¬ 
namical system given by Eq. (|T|) (see 11811 for a detailed de¬ 
scription of the experimental realization of the Rossler-like 
oscillator, and [3Cj, 3l]] for previous realizations in network 
configurations). As in the evaluation of the MSFs illustrated 
in Fig. [3 we consider the system coupled through the x and 
then through the y variable, the cases being paradigmatic ex¬ 
amples of class III and class II systems, respectively. 


The experimental setup shown in Fig. [3 consists of an 
electronic array (EA), a multifunction data card (DAQ), and 
a personal computer (PC). The EA comprises 6 Rossler-like 
electronic circuits forming the spiderweb network with one 
central node and 5 peripheral nodes. Each node has an indi¬ 
vidual electronic coupler controlled by a digital potentiometer 
(XDCP), which is adjusted by a digital output signal (DO) 
coming from ports PO.O and P0.1. Port PO.O is used to set the 
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value of the coupling resistance (adequately scaled to corre¬ 
spond to the values of a) and PO. 1 increases or decreases the 
value of the resistance through a voltage divisor (the resolu¬ 
tion allowing for 100 discretized steps). The full experimen¬ 
tal process is controlled with a virtual interface developed in 
Labview 8.5, that can be considered as a state machine. The 
experimental procedure is realized as follows. First, a is set to 
zero, after a waiting time of 500 ms (roughly corresponding 
to 600 cycles of the autonomous systems), the output signals 
from the 6 circuits are acquired by the analog ports (AI 0; 
AI 1; ...; AI 5). Once the dynamics of the whole ensemble 
is recorded, the value of a is increased by one step, and the 
signals are again stored in the PC for further analysis. This 
process is repeated until the maximum value of a is reached. 

Figure [7] shows the synchronization error of the whole net¬ 
work calculated as N .^_^ \ x i ~ x j\> w h ere the nor¬ 

malizing factor corresponds to the total number of oscilla¬ 
tor pairs in the network. We repeat the experiment from 10 
different initial conditions and compute the average synchro¬ 
nization error across the realizations. Whatever the initial 
condition, the 6-oscillator network exhibits strongly intermit¬ 
tent dynamics which manifests itself by intermittent switches 
between S and L in every oscillator node. Therefore, we 
compute the synchronization error in three different ways: I) 
the total error of the whole time series regardless of visited 
states, T, blue (dark gray) dashed line II) the error during the 
time windows where all oscillators exhibit S dynamics, red 
(gray) dash-dotted line, and III) the error during time win¬ 
dows, where all oscillators represent L dynamics, green (light 
gray) continuous line. The synchronization errors for these 
three cases when the coupling is introduced through the x 
variable is shown in Fig. [7] (a). In the topology under study, 
the ratio of the largest to the smallest nonzero eigenvalue of 
the Laplacian matrix is Xn /^2 = 2.519, that is smaller than 
the ratio of the largest to the smallest zeros of the MSF of 
L vi/v\ = 3.97 (and obviously smaller than the conspicu¬ 
ously larger analogous ratio for S). There is thus a range of 



FIG. 6. (Color online) Experimental setup. On the left, schematic 
representation of the coupling topology of the 6-circuit network. The 
coupling is adjusted by a digital potentiometers X9C104, whose pa¬ 
rameters C u /d (Up/Down resistance) and C step (increment of the re¬ 
sistance at each step) are controlled by a digital signal coming from 
a DAQ Card, P0.0-P0.1 respectively. The outputs of the circuit are 
sent to a set of voltage followers that act as a buffer and, then, sent to 
the analog ports (AI 0 ; AI 1; ... ; AI 5) of the same DAQ Card. The 
whole experiment is controlled by a PC using a Labview Software. 


values of the coupling strength a for which all eigenmodes 
are predicted to be in the stability region (i.e., 0 -A 2 > 
and erAjv < ^2 hold simultaneously), and the system is syn- 
chronizable. Indeed, we observe that, when er increases, com¬ 
plete synchronization of the whole network is achieved around 
a = 0.165 132)]. The experimental system becomes unstable 
(the oscillations suddenly disappear and the system reaches 
a fixed point) above a ~ 0.25, so our results are most rele¬ 
vant in relation to the eigenmode corresponding to A 2 and the 
crossing of the first MSF zero (given the Xn /^2 ratio, at the 
time this eigenmode crosses only slightly the stability region, 
the other eigenmodes are guaranteed to be stable). Figure [7] 
(b) shows the same results for the case of y-coupling. In this 
case the system is indeed class II, so again we focus on the 
crossing of the eigenmode corresponding to A 2 beyond the 
(only) zero of the MSF. The network reaches synchronization 
for a > 0.0435 and does not leave the synchronized manifold 
for larger values of the coupling strength, as predicted by the 
MSF. 

At this point, it is interesting to analyze the influence of 
the parameter mismatch on the emergence of synchroniza- 


(a) 



a 


(b) 
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FIG. 7. (Color online) Experimental results: synchronization er¬ 
ror as a function of the coupling strength, (a) Synchronization 
error (see main text for a definition) for a network of oscillators cou¬ 
pled through the x variable (class III system) for whole time series T 
(blue, dark gray, middle curves), time windows ss where all oscilla¬ 
tors exhibit S dynamics (red, gray, lower curves), and time windows 
ll where all oscillators have L dynamics (green, light gray, upper 
curves). Insets: the curves zoomed around the smallest synchroniza¬ 
tion error achieved in the system, the results obtained from all 10 
individual realizations (dotted lines) superimposed on the averages 
(continuous lines), (b) Analogous curves for a network of oscillators 
coupled through the y variable (class II system). 
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tion. First, we study the effect of parameter uncertainties in 
the MSFs. All parameters in our system (aq, a 2 , < 23 , /3, T, 7 , 
S and /i in Eq. (Q3) are resistances, capacitances, and products 
thereof (see II1811 for the details). The insets of Fig. [5] show 
the effect of the uncertainties in these parameters (taking the 
resistor tolerance to be 1% and the capacitor tolerance to be 
10%) on 50 realizations for each possible case (x-coupling 
for S and L in panel (a), and //-coupling for S and L in panel 
(b)). We zoom around the first zero of each MSF, as this is the 
most relevant region for the comparison with experiment. De¬ 
viations from the noiseless case of the order of 20% or higher 
are frequently seen. Nevertheless, the MSFs for S and L do 
not overlap despite the parameter fluctuations and therefore 
are easily distinguishable. This is in a good qualitative agree¬ 
ment with our experiments. The insets of Figs. [7](a)-(b) show 
the results zoomed around the transition to the synchronized 
regime, indicating the values of the 10 individual realizations 
and their corresponding average, showing a relatively small 
variability that qualitatively agrees with our numerical predic¬ 
tions on the effects of a parameter mismatch as shown in Fig. 
0 

Finally, Fig. [8] shows the fraction of time that the intermit¬ 
tent / system spends in ss, red (gray) triangles, ll, green (light 
gray) squares, and si, black (dark gray) diamonds, time win¬ 
dows. One can see the qualitative agreement with the case of 
two coupled oscillators shown in Fig. Q](a). 


(a) 



(b) 



FIG. 8. (Color online) Experimental results: fraction of time that 
intermittent system spends in ss, red (gray) triangles, ll, green 
(light gray) squares, and si, black (dark gray) diamonds, time 
windows as a function of a. (a) Network coupled through the x 
variable and (b) network coupled through the y variable. The results 
confirm qualitatively the numerically predicted behavior shown in 
Fig.il 


CONCLUSIONS 

We have shown how the MSF approach can be used for 
the analysis of synchronization in ensembles of multistable 
chaotic systems. To introduce our methodology, we have used 
a Rossler-like oscillator with two coexisting chaotic attractors. 
The existence of interwoven basins of attraction hinders the 
prediction of the system asymptotic behavior in the presence 
of noise or parameter uncertainties. Under this framework, 
we have analyzed the synchronization regimes of two bidi¬ 
rectionally coupled chaotic oscillators and showed the enor¬ 
mous complexity of the basins of attraction leading to diverse 
dynamics of the whole system, as well as the fact that the 
coexisting attractors became synchronized at different cou¬ 
pling strengths, manifesting various synchronization types for 
a given coupling strength. After that, we have proposed the 
use of the MSF arguments for multistable systems, which pro¬ 
vided information on synchronizability of a given network 
of multistable oscillators. Specifically, the MSF shows un¬ 
der what conditions a network of multistable systems can 
synchronize for a given range of topology spectra and cou¬ 
pling strengths, whatever might be the attractor dynamics to 
which different oscillators become locked in the presence of 
intermittency. Even though in this work we have focused on 
chaotic dynamics, it is important to point out that our ap¬ 
proach, as based on a MSF reasoning, can also be used to 
study synchronization of multistable periodic systems (a MSF 
approach to monostable periodic, as well as chaotic, dynamics 
appears for example in |33|]). Finally, we have experimentally 
demonstrated the feasibility of the MSF approach with a net¬ 
work of oscillating circuits in a heavily intermittent regime, 
and shown that the predictions, as in the case of monostable 
systems, affected by small uncertainties, were nonetheless 
very useful from a qualitatively point of view, showing the ro¬ 
bustness of the proposed methodology to noise and parameter 
mismatch. 

The proposed methodology is of a special interest for a 
stability analysis of synchronization in multistable systems 
during intermittency, since any monostable approach to syn¬ 
chronization is bound to fail in that regime, and knowing the 
coupling strength or the topological modifications required to 
maintain complete synchronization under any possible attrac¬ 
tor dynamics is especially useful in that scenario. The MSF 
arguments in the context of multistable systems provide a 
generic tool to understand complete synchronization of real 
multistable systems such as those occurring in laser physics 
0, genetics (§] or cell signaling 
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APPENDIX - COMPUTATION OF LYAPUNOV EXPONENTS 

The computation of Lyapunov exponents of coupled 
Rossler-like systems throughout this paper has been per¬ 
formed using the well established method first proposed and 
justified theoretically in Refs. 0. The reader may, how¬ 
ever, raise the objection that the equations of motion of the 
system under study are given by a non-smooth vector field, 
whose first derivate presents a discontinuity, and therefore the 
validity of our results is not clear at this point. Indeed, the 
equation of motion for the z variable in Eq. O contains a 
piecewise linear function g{x), whose derivative g'{x) is such 
that linLj._,. 3 - g'(x) ^ lim x _ >3 + g'(x). Even if the mathemat¬ 
ical problem of the existence of Lyapunov exponents for a 
sufficiently general dissipative dynamical system remains un¬ 
solved to this very day |35i| . from a more practical perspective, 
the difficulty of computing quantities that are based on the lin¬ 
earized variational equations of a dynamical system whose Ja¬ 
cobian matrix is undefined in the phase-space plane x = 3 is 
already apparent. So there are grounds for this objection, even 
if in practice the Lyapunov spectrum of our systems converge 
to well defined asymptotic values. On the other hand, a more 
practically minded reader would point out that the plane x = 3 
is a zero Lebesgue measure set, and that even in the discretized 
world of numerical analysis the chances that we get states cod¬ 
ified with double precision such that x = 3.000000... are 
very small. The argument is that if indeed such states ever 
occur, their contribution to the Lyapunov exponents (which 
are effectively computed as long time averages across phase- 
space orbits) will be negligible. We will see that this latter 
position turns out to be vindicated by a close examination of 
the issue. 

To avoid the above mentioned non-smoothness, the idea of 
using a slightly modified dynamical system, where the dis¬ 
continuity at x = 3 is bypassed at sufficiently small scales in 
such a way that it does not affect too much the estimation of 
the Lyapunov spectrum, is one that naturally comes to mind. 
As a matter of fact, our choice of a dynamical system was ulti¬ 
mately determined by the fact that the piecewise Rossler sys¬ 
tem is known for its robustness and experimental accessibility 
when implemented in electronic circuits. There, the disconti¬ 
nuity in the slope of g(x) is obtained by using a semiconductor 
diode placed in a voltage divider, whose I — V characteristics 
is supposed to be such that I > 0 for V = Vd and 7 = 0 for 
V <Vd- But of course, this is known to be a crude simplifica¬ 
tion of the diode behavior, as statistical physics models con¬ 
sidering the transport of charge carriers across the depletion 
layer at the semiconductor p — n junction result in a smooth 
function in the macroscopic I — V behavior. One possibility 
would be to include equations based on the Shockley diode 
law or more realistic mathematical models of a diode, and ob¬ 
tain the relevant g(x) of what would be a considerably more 
complicated smooth dynamical system. Instead, we want to 


explore the possibility of simply replacing g(x) with a per¬ 
turbed smooth function h(x) defined as follows 

( 0 x < 3, 

h(x) = < p(x) 3 < x < 3 + 5x, (4) 

[ p(x — 3) x>3 + Sx, 

where Sx <C 3. Here, p(x) can be a polynomial of a degree so 
high to have as many continuous derivatives as one wishes at 
x = 3 and x = 3 + Sx, or a more sophisticated interpolating 
function. For the sake of simplicity, let us take p(x) to be a 
third degree polynomial. Four undetermined constants guar¬ 
antee that we can satisfy four minimal conditions to provide 
continuity of h(x) and h'(x) at both x = 3 and x = 3 + Sx. 

The resulting polynomial p(x) = — {x — 3) 3 + (x — 3) 2 

is shown in Fig. [9] for Sx = 0.010 (blue line), 0.005 (red 
line) and 0.002 (green line). We also include the limiting case 
Sx = 0 (black line), for which h(x) = g(x). 

Next, we check how the estimated Lyapunov exponents 
of the Rossler oscillator are affected by replacing g(x) with 
h[x) in the equations of motion. We expect the effect to be 
smaller and smaller as Sx decreases. This expectation is in¬ 
deed matched by the results, as can be seen in the upper panel 
of Fig. [TO] Here, the largest Lyapunov exponents correspond¬ 
ing to attractors S (magenta diamonds) and L (gray circles) 
are shown. They result from the integration of the system 
across a time window comprising roughly 30000 cycles. For 
large Sx, the change in the dynamics significantly affects the 
Lyapunov exponents, as it is seen quite clearly in the fact that 
S for some choices of Sx becomes a regular (i.e. non-chaotic) 
attractor with a zero maximum Lyapunov exponent. On the 
other hand, for any choice of Sx such that Sx < 0.01 all max¬ 
imum Lyapunov exponent estimates converge, within our ac¬ 
curacy, to the same values, and that is also the case for the rest 
of the Lyapunov spectrum (not shown here). 

The previous results can be simply related to the fact that 



FIG. 9. (Color online) Function g(x) around x = 3 (black line, 
Sx = 0), and function h(x) for three different Sx values (blue 
line, Sx = 0.010; red line, Sx = 0.005; green line, Sx = 0.002). 
For sufficiently small Sx, h(x) can be made as close as needed to 
g(x) without discontinuities in the first derivative. Higher degree 
polynomials could be used in a completely analogous manner to ap¬ 
proximate g(x) with a higher number of continuous derivatives at 
x = 3 and x = 3 + Sx. 
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for small Sx a small fraction of the phase space points fall in¬ 
side the region where h(x) differs from g(x). Indeed, we see 
in Fig. [TO] (b) that this fraction decreases as Sx is reduced as 
a power law with a characteristic exponent close to 1. There¬ 
fore, replacing g[x) with h(x) with a sufficiently small Sx 
solves the issue of the discontinuity in the Jacobian of the sys¬ 
tem giving Lyapunov exponent estimates independent of the 
precise values of Sx. 

Finally, we consider the possibility of taking the limiting 
case Sx = 0 ( h(x ) = g(x)), which is tantamount to simply 
ignoring the discontinuity in the Jacobian. In practice, we use 
the criterion that if a phase space point happens to be such 
that x = 3, we consider the derivative of the z component of 
the vector field with respect to x to be its left derivative (i.e. 
the derivative of g(x) is set to lim 2 ._ ) , 3 - g'(x) = 0), which 
is also our choice throughout the paper. The corresponding 
Lyapunov exponent estimates are shown in Fig. [TO] (a) as the 
dashed lines (magenta for S, gray for L). They coincide with 
the values obtained by replacing g(x) with h(x) for a suffi¬ 
ciently small Sx. We could also choose to assign as derivative 
of g{x) at x = 3 half of the time the left derivative, half of 
the time the right derivative, or any other reasonable criterion. 
Ultimately, this choice is inconsequential, as the fraction of 
points that with x = 3 is zero for both L and S, none of the 


(a) 




FIG. 10. (Color online) Maximum Lyapunov exponents as a 
function of 8x, and fraction of phase space points for which 

x £ [3, 3 + 5x\. In (a), the largest Lyapunov exponents corre¬ 
sponding to attractors S (magenta diamonds) and L (grey circles) 
as functions of Sx. The dashed lines in the background correspond 
to the estimates obtained with Sx = 0 ( h(x ) = g{x)). In (b), the 
fraction of phase space points within the interval [3, 3 + <5ir], where 
h(x) + g(x). 


considered 30 million phase-space points has an x coordinate 
that is exactly 3 (in the double precision floating point for¬ 
mat). As expected from the results shown in Fig. |T0](b), this 
will also be the case for the slightly perturbed system with 
h(x ) in lieu of g{x) if Sx is so small that the fraction of points 
with x £ [3, 3 + Sx] times the total number of orbit points 
considered is considerably smaller than a number on the order 
of unity. Even if a few points in hundreds of thousands or mil¬ 
lions of phase-space points were affected by the discontinuity 
in the Jacobian or by the small modification in the dynamics 
introduced by the polynomial p(x), the effect on long time 
averages along phase-space orbits, such as those Lyapunov 
exponent estimates are based upon, would be negligible. 

In conclusion, in order to avoid the discontinuity in the Ja¬ 
cobian of our dynamical system, we can compute the Lya¬ 
punov exponents using a slightly modified dynamical system. 
Doing this, we must ensure that the modification should be 
sufficiently small so that the results become independent of 
the size of the phase-space region whose dynamics is modified 
and of the modification form (independent of the size of the 
modified region Sx and of the degree of the polynomial p(x) 
in the example above). But then, the resulting Lyapunov expo¬ 
nent estimates coincide with those obtained from the original 
dynamical system. 
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